rm(list = ls())

library(ctgt)
library(hommel)
library(faux)

## metabolites
n <- 50
m <- 100
rho = 0
metab <- paste("x", 1:m, sep = "")
### pathway list
pathEle <- list()
for(i in 1:10){
  pathEle[[i]] <- paste("x", (10*(i-1)+1):(10*i), sep = "")
}


pathList = list(manystrong = unlist(pathEle[1]),
                fewstrong = c(metab[11:15], unlist(pathEle[9:10]) ) ,
                moreweak = unlist(pathEle[3:8]),
                manyweak = unlist(pathEle[3:5]),
                fewweak = c(metab[16:20], unlist(pathEle[9:10]) ),
                null = unlist(pathEle[9]) )


lp = length(pathList)
run = 1000
res_st_mn9_100 = matrix(0, lp, run)
res_gt_mn9_100 = matrix(0, lp, run)

### data set
for(s in 1:run){
  set.seed(1234+s)
  #  X = matrix(0,n,m)
  #  colnames(X) <- metab
  #  for(i in 1:n){
  #    X[i,] <- arima.sim(list(order = c(1,0,0), ar = 0.9), n = m)
  #  }
  
  X = as.matrix(rnorm_multi(n, m, 0, 1, rho, varnames = metab))
  y <- rbinom(n = n, size = 1, prob = 0.5)
  
  X[y == 1, 1:10] <- X[y == 1, 1:10] + 0.95
  X[y == 1, 11:15] <- X[y == 1, 11:15] + 1.05
  
  X[y == 1, 21:80] <- X[y == 1, 21:80] + 0.55
  X[y == 1, 16:20] <- X[y == 1, 16:20] + 0.7
  
  
  
  pvs = sapply(metab, function(j) gt2(y,X,hyps = j))[1,]
  hom = hommel(pvs, simes = TRUE)
  
  for(pathid in 1:length(pathList)){
    lt = localtest(hom, pathList[[pathid]], tdp=0)
    res_st_mn9_100[pathid, s] = ifelse(lt<0.05, "reject", "not reject")
    res_gt_mn9_100[pathid, s] = actgt(y,X,metab,hyps = pathList[[pathid]], maxit = 5000, alpha = 0.05 )[1]
  }
  
}

save(res_st_mn9_100, file = "st_mn0_100.RData")
save(res_gt_mn9_100, file = "gt_mn0_100.RData")



rm(list = ls())

library(ctgt)
library(hommel)
library(faux)

## metabolites
n <- 50
m <- 100
rho = 0.9
metab <- paste("x", 1:m, sep = "")
### pathway list
pathEle <- list()
for(i in 1:10){
  pathEle[[i]] <- paste("x", (10*(i-1)+1):(10*i), sep = "")
}


pathList = list(manystrong = unlist(pathEle[1]),
                fewstrong = c(metab[11:15], unlist(pathEle[9:10]) ) ,
                moreweak = unlist(pathEle[3:8]),
                manyweak = unlist(pathEle[3:5]),
                fewweak = c(metab[16:20], unlist(pathEle[9:10]) ),
                null = unlist(pathEle[9]) )


lp = length(pathList)
run = 1000
res_st_mn9_100 = matrix(0, lp, run)
res_gt_mn9_100 = matrix(0, lp, run)

### data set
for(s in 1:run){
  set.seed(1234+s)
  #  X = matrix(0,n,m)
  #  colnames(X) <- metab
  #  for(i in 1:n){
  #    X[i,] <- arima.sim(list(order = c(1,0,0), ar = 0.9), n = m)
  #  }
  
  X = as.matrix(rnorm_multi(n, m, 0, 1, rho, varnames = metab))
  y <- rbinom(n = n, size = 1, prob = 0.5)
  
  X[y == 1, 1:10] <- X[y == 1, 1:10] + 0.95
  X[y == 1, 11:15] <- X[y == 1, 11:15] + 1.05
  
  X[y == 1, 21:80] <- X[y == 1, 21:80] + 0.55
  X[y == 1, 16:20] <- X[y == 1, 16:20] + 0.7
  
  
  
  pvs = sapply(metab, function(j) gt2(y,X,hyps = j))[1,]
  hom = hommel(pvs, simes = TRUE)
  
  for(pathid in 1:length(pathList)){
    lt = localtest(hom, pathList[[pathid]], tdp=0)
    res_st_mn9_100[pathid, s] = ifelse(lt<0.05, "reject", "not reject")
    res_gt_mn9_100[pathid, s] = actgt(y,X,metab,hyps = pathList[[pathid]], maxit = 5000, alpha = 0.05 )[1]
  }
  
}

save(res_st_mn9_100, file = "st_mn9_100.RData")
save(res_gt_mn9_100, file = "gt_mn9_100.RData")



rm(list = ls())

library(ctgt)
library(hommel)
library(faux)

## metabolites
n <- 50
m <- 150
rho = 0
metab <- paste("x", 1:m, sep = "")
### pathway list
pathEle <- list()
for(i in 1:10){
  pathEle[[i]] <- paste("x", (10*(i-1)+1):(10*i), sep = "")
}


pathList = list(manystrong = unlist(pathEle[1]),
                fewstrong = c(metab[11:15], unlist(pathEle[9:10]) ) ,
                moreweak = unlist(pathEle[3:8]),
                manyweak = unlist(pathEle[3:5]),
                fewweak = c(metab[16:20], unlist(pathEle[9:10]) ),
                null = unlist(pathEle[9]) )


lp = length(pathList)
run = 1000
res_st_mn9_100 = matrix(0, lp, run)
res_gt_mn9_100 = matrix(0, lp, run)

### data set
for(s in 1:run){
  set.seed(1234+s)
  #  X = matrix(0,n,m)
  #  colnames(X) <- metab
  #  for(i in 1:n){
  #    X[i,] <- arima.sim(list(order = c(1,0,0), ar = 0.9), n = m)
  #  }
  
  X = as.matrix(rnorm_multi(n, m, 0, 1, rho, varnames = metab))
  y <- rbinom(n = n, size = 1, prob = 0.5)
  
  X[y == 1, 1:10] <- X[y == 1, 1:10] + 0.95
  X[y == 1, 11:15] <- X[y == 1, 11:15] + 1.05
  
  X[y == 1, 21:80] <- X[y == 1, 21:80] + 0.55
  X[y == 1, 16:20] <- X[y == 1, 16:20] + 0.7
  
  
  
  pvs = sapply(metab, function(j) gt2(y,X,hyps = j))[1,]
  hom = hommel(pvs, simes = TRUE)
  
  for(pathid in 1:length(pathList)){
    lt = localtest(hom, pathList[[pathid]], tdp=0)
    res_st_mn9_100[pathid, s] = ifelse(lt<0.05, "reject", "not reject")
    res_gt_mn9_100[pathid, s] = actgt(y,X,metab,hyps = pathList[[pathid]], maxit = 5000, alpha = 0.05 )[1]
  }
  
}

save(res_st_mn9_100, file = "st_mn0_150.RData")
save(res_gt_mn9_100, file = "gt_mn0_150.RData")











rm(list = ls())

library(ctgt)
library(hommel)
library(faux)

## metabolites
n <- 50
m <- 150
rho = 0.9
metab <- paste("x", 1:m, sep = "")
### pathway list
pathEle <- list()
for(i in 1:10){
  pathEle[[i]] <- paste("x", (10*(i-1)+1):(10*i), sep = "")
}


pathList = list(manystrong = unlist(pathEle[1]),
                fewstrong = c(metab[11:15], unlist(pathEle[9:10]) ) ,
                moreweak = unlist(pathEle[3:8]),
                manyweak = unlist(pathEle[3:5]),
                fewweak = c(metab[16:20], unlist(pathEle[9:10]) ),
                null = unlist(pathEle[9]) )


lp = length(pathList)
run = 1000
res_st_mn9_100 = matrix(0, lp, run)
res_gt_mn9_100 = matrix(0, lp, run)

### data set
for(s in 1:run){
  set.seed(1234+s)
  #  X = matrix(0,n,m)
  #  colnames(X) <- metab
  #  for(i in 1:n){
  #    X[i,] <- arima.sim(list(order = c(1,0,0), ar = 0.9), n = m)
  #  }
  
  X = as.matrix(rnorm_multi(n, m, 0, 1, rho, varnames = metab))
  y <- rbinom(n = n, size = 1, prob = 0.5)
  
  X[y == 1, 1:10] <- X[y == 1, 1:10] + 0.95
  X[y == 1, 11:15] <- X[y == 1, 11:15] + 1.05
  
  X[y == 1, 21:80] <- X[y == 1, 21:80] + 0.55
  X[y == 1, 16:20] <- X[y == 1, 16:20] + 0.7
  
  
  
  pvs = sapply(metab, function(j) gt2(y,X,hyps = j))[1,]
  hom = hommel(pvs, simes = TRUE)
  
  for(pathid in 1:length(pathList)){
    lt = localtest(hom, pathList[[pathid]], tdp=0)
    res_st_mn9_100[pathid, s] = ifelse(lt<0.05, "reject", "not reject")
    res_gt_mn9_100[pathid, s] = actgt(y,X,metab,hyps = pathList[[pathid]], maxit = 5000, alpha = 0.05 )[1]
  }
  
}

save(res_st_mn9_100, file = "st_mn9_150.RData")
save(res_gt_mn9_100, file = "gt_mn9_150.RData")
